Sparse regression problems arise often in various applications, e.g., model selection, compressive sensing, EEG source localisation and gene modelling [1], [2]. One of the Bayesian approaches to force the coefficients being zeros is the spike and slab prior [3]: each component is modelled as a mixture of spike, that is the delta-function in zero, and slab, that is some vague distribution. Following the Bayesian approach, latent variables that are indicators of spikes are added to the model [4] and the relevant distribution is placed over them [5]. In this model each component is modelled to be spike or slab independently. However, in many applications nonzero elements tend to appear in groups forming an unknown structure: wavelet coefficients of images are usually organised in trees [6], chromosomes have a spatial structure along the genome [2]. We propose an extension of the spike and slab model by imposing a hierarchical Gaussian process (GP) prior on the latent variables. Such hierarchical prior allows to model spatial structural dependencies for coefficients that can evolve in time. The new model is flexible as spatial and temporal dependencies are decoupled by different levels of the hierarchical GP prior.